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We investigate magnetic properties of strongly interacting four component spin-3/2 ultracold 
fermionic atoms in the Mott insulator limit with one particle per site in an optical lattice with 
honeycomb symmetry. In this limit, atomic tunneling is virtual, and only the atomic spins can 
exchange. We find a competition between symmetry breaking and liquid like disordered phases. 
Particularly interesting are valence bond states with bond centered magnetizations, situated between 
the ferromagnetic and conventional valence bond phases. In the framework of a mean-field theory, 
we calculate the phase diagram and identify an experimentally relevant parameter region where a 
homogeneous SU(4) symmetric AfHeck-Kennedy-Lieb-Tasaki-like valence bond state is present. 


I. INTRODUCTION 


Mott insulators with antiferromagnetic spin correla¬ 
tions are in the center of interest in condensed matter 
physics due to their relation to high-Tc superconductiv¬ 
ity [1-3]. They are also intensively investigated in quan¬ 
tum information, because of their potential applicability 
for quantum computing [4]. In particular, these systems 
can be used to realize states required for measurement- 
based quantum computation (MBQC) [5-9]. In addi¬ 
tion they may exhibit nontrivial topology [10-13]; ac¬ 
cordingly, they play an important role in the studies of 
the topological states of matter. For experimental studies 
of these phenomena ultracold atomic systems provide one 
of the most promising and efficient playgrounds. There is 
a still ongoing progress in ground breaking experiments 
with ultracold atoms in order to realize quantum emu¬ 
lators of magnetic systems [14]. The main advantage is 
that in these systems there is an unprecedented control 
over the parameters describing almost every feature of 
their physics [15-21]. The atoms are trapped optically; 
both the potential height and the lattice periodicity can 
be adjusted by tuning the amplitude, phase and wave¬ 
length of the lasers. Even the geometry of the lattice 
can be changed in situ [22]. Further advantage of such 
systems is that interaction between the atoms can be 
controlled in a wide range through the access of various 
scattering resonances [23-26]. As a result of this versa¬ 
tility, many antiferromagnets, either encountered in real 
materials, or proposed by theorists for academic inter¬ 
est, can now be potentially realized [27-37]. In the first 
experiments the Mott insulator state was realized with 
a dilute gas sample of alkaline atoms loaded to an op¬ 
tical lattice formed by counterpropagating laser beams 
[38-40]. Later, trapping of higher spin alkalies [41, 42] 
and cooling of alkaline-earth-metal atoms to quantum de¬ 
generacy [43, 44] has opened the way to Mott insulators 


with higher spin atoms [45, 46]. Most recently, the first 
steps were made towards the study of the direct effects 
of spin-exchange interactions [47-49] 

The word ’antiferromagnet’ refers to a state of mat¬ 
ter, where the total magnetization of the sample is zero, 
but magnetic correlations differ from that of the triv¬ 
ial paramagnetic phase. The simplest possible antiferro¬ 
magnetic state has Neel order: in a square lattice in 2D, 
for instance, opposing spins are arranged in a checker¬ 
board configuration. Such a state is symmetry breaking, 
and an order parameter can be introduced, which is the 
magnetization of the sublattice formed by every second 
site. According to the Mermin-Wagner-Hohenberg the¬ 
orem [50, 51] symmetry breaking cannot take place in 
two dimensions at finite temperature, while in one di¬ 
mension fluctuations destroy long-range order even at 
T = 0 [52]. Therefore, another completely different an¬ 
tiferromagnetic state can be expected in low dimensions: 
the singlet covering of the lattice, where pairs of sites 
form a two-particle spin-singlet state [53]. Similar state 
was discovered by Majumdar and Ghosh in a ID model 
[54], where the ground state is a periodic lattice of in¬ 
dependent singlet pairs. A translational invariant gen¬ 
eralization of such a valence-bond-solid (VBS) state was 
introduced by Anderson [2] to describe frustrated anti¬ 
ferromagnetism on a triangular lattice; it is called the 
resonating valence bond (RVB) state. 

The direct consequences of an antiferromagnetic spin 
exchange can be studied most clearly in the Mott insula¬ 
tor state of matter [1, 53, 55-57]. In this case the charge 
degree of freedom is frozen out — e.g. as a consequence 
of a strong repulsive interaction — and the low energy 
physics of the system is governed by the spin degrees 
of freedom, which in turn are governed by an effective 
Heisenberg-like exchange interactions. However, as men¬ 
tioned above, it is not obvious, particularly in 2D, that 
even in the case of strong antiferromagnetic exchange in¬ 
teractions, a Neel-type spin ordered state will emerge. 
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Instead, valence bonds can form and the SU(2) spin ro¬ 
tational invariance can remain preserved. Nowadays the 
concept of the valence bond picture is the most com¬ 
monly used to describe the underlying order of the vari¬ 
ous spin liquid states. With the help of the valence bond 
picture a series of phenomena with nontrivial magnetic 
origin can be described like various exotic VBS states 
[27, 29, 30, 58-60], or even chiral spin liquid (CSL) states 
with non-trivial topology [13, 27, 30, 31, 59]. In Ref. [61] 
it was shown that a topological charge-Haldane state can 
emerge on a spin-3/2 fermionic ladder. The SU(4) sym¬ 
metric spin-3/2 system has been studied intensively in 
the last few years [29, 58, 62-65], mostly on the square 
lattice. Although some recent numerical results suggest 
that a weak SU(4) symmetry breaking order can emerge 
[62], it is mostly accepted that the ground state is a 
VBS state consisting disconnected RVB plaquettes. In 
Ref. [63] a pure SU(4) Heisenberg model was studied on a 
honeycomb lattice; the authors found that a spin-orbital 
liquid phase can emerge in the system which collapse 
into a tetramerized VBS-like state in the presence of next 
nearest neighbor exchange [66]. 

The states, which can be suitable as universal resource 
for MBQC, could be searched among these spin liquid¬ 
like states, which preserve the usual SU(2) spin rotational 
invariance. The classification of these states, just as the 
quest for physical systems realizing them and harvesting 
their potential use as universal resource states, is a big 
challenge of quantum information theory. A fundamental 
requirement is that they have to be the unique ground 
states of a Hamiltonian with gapped spectrum and short- 
range interactions in order to assure their robustness [67]. 
One of the most promising states is the AKLT-state in¬ 
troduced by Affleck, Kennedy, Lieb and Tasaki who pro¬ 
posed several models in ’dimension one and more’ where 
the ground state is a unique VBS [68, 69]. In the contem¬ 
porary technical ” slang” these are states that are exactly 
described by the Matrix Product States (MPS) or Pair 
Entangled Projected States (PEPs) with the lowest pos¬ 
sible (nontrivial) local dimension (cf. [14]). 

However, the ’parent’ Hamiltonians of these ideal¬ 
ized states are not easy to realize as low energy ef¬ 
fective Hamiltonian of a fermionic or bosonic Mott in¬ 
sulator with SU (2) invariant interactions. Let us re¬ 
mind the reader that when the internal states of par¬ 
ticles correspond to the degenerate manifold of the hy- 
perfine atomic level A, i.e. IE, mi?) of an ultracold 
spinor gas, the resulting Hamiltonians, in the absence 
of symmetry breaking fields, must be SU{2) symmetric. 
The corresponding spin Hamiltonians are given then by 
powers of the nearest neighbor Heisenberg interactions 
H — j> where Uk are numbers that 

depend on scattering lengths (for a review see for instance 
[70] for E = 1,3/2, 2, 5/2, and for general Fermi systems 
c.f. [27, 71, 72]). 

The situation for the Mott insulators with one particle 
per site might be summarized as described below; the 
case of Mott insulators with two or more particles per 


site is much more complex (cf. [70]). 

• The original AKLT state was proposed for a spin 

E = 1 in ID. E = 1 can s-wave collide in the 
Ftot = 0, 2 channels, and thus are characterized in 
general by the two distinct scattering lengths. The 
effective Hamiltonian in the super-exchange limit 
reads then H = J2kZo ^ki^i ■ or equiva¬ 

lently H = E(z,i> Es=o, 2 fe'^s(S* -f Sj)], where 
qs are numbers, and Ps are projections on the 
total spin S. If one could control the two scat¬ 
tering lengths independently and arbitrarily, one 
could realize the case when H = J2(i j) ’^ 2 (Si+Sj)], 
which in ID corresponds exactly to the AKLT case. 
Unfortunately, such control of scattering lengths 
nowadays is hardly possible - see ref. [70] for de¬ 
tails. 

• The Hamiltonian H = 

E = 1 particles is known as biquadratic-bilinear 
Hamiltonian and at least in ID has been studied 
very intensively (cf. [73-77] and references therein). 
It is known that the antiferromagnetic regime ex¬ 
hibits a robust gapped phase that does not break 
the SU{2), the celebrated Haldane phase [78, 79]. 

• For E = 2 bosons there are three possible s-wave 

scattering channels and three scattering length, re¬ 
spectively. The effective Hamiltonian has the form 
H = Es=o.2,4[9siPs(Si -f Sj)], and can in 

principle be reduced to H = J2{i j) ^ 4 (Si + Sj)] by 
adjusting the scattering lengths. That would corre¬ 
spond to AKLT model on the square lattice or the 
3D lattice with coordination number 4. Again, in 
practice the necessary control of scattering length 
is not possible. 

• For E = 3/2 fermions there are two possible s-wave 

scattering channels and, consequently, two scatter¬ 
ing length, respectively. The effective Hamiltonian 
has the form H = Es=o.2bs7^s(S* -f S^-)], 

and cannot be by any means reduced to H = 

+ Sj)], which on the honeycomb lat¬ 
tice in 2D corresponds to the AKLT model. Similar 
situation holds for higher half-integer spin. 

• The effective model with E = 3/2 was studied ex¬ 
tensively, and already few years ago there has been 
a lot of progress in understanding the special prop¬ 
erties of E = 3/2 and E = 5/2 Fermi gases. In 
spin-3/2 systems with contact interaction, Wu et 
al. realized that a generic SO(5) symmetry exists 
[80] . These authors also found novel competing or¬ 
ders [81, 82], suggesting a quartetting ph.a,se and the 
s -wave quintet Cooper pairing phase. 

Let us recapitulate: the parent Hamiltonian of the 
two-dimensional generalization of the AKLT-state is a 
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spin-Hamiltonian of spin-3/2 operators on a two- Hubbard model, 
dimensional honeycomb lattice: 


= 53^^3(S. + S,), 

(iJ) 


( 1 ) 




- V r r 

i,o',0 


( 2 ) 


where 7^3(Si -1- Sj) projects to the total spin-3 subspace: 
(Si -1- S_,)^ = 3(3-1- 1). In the next section we will demon¬ 
strate in detail that the effective spin Hamiltonian de¬ 
rived from a usual Hubbard-like model does not contain 
the V 3 projector. Although, in a more general model this 
term would appear, its role remains always secondary, 
unless we assume ’’all mighty” control over the scatter¬ 
ing length - typically this term is generated by weaker 
interactions: as a perturbation it can change the ground 
state even drastically, but it will never become dominant. 

In this paper we show that a spin-3/2 ultracold 
fermionic system loaded into a two-dimensional hon¬ 
eycomb lattice has a ground state similar to the two- 
dimensional generalization of the AKLT-state for an ex¬ 
tended parameter range of the coupling constants de¬ 
scribing the on-site fermion-fermion interaction. Our 
analysis is based on a mean-field study of the system with 
a suitable ansatz to describe the coexistence of or even 
the competition between site- and bond-centered spin or¬ 
ders. We find that the AKLT-like homogeneous state 
is the lowest energy solution in an extended experimen¬ 
tally reachable parameter region. This state competes on 
one hand with the usual site ordered homogeneous fer¬ 
romagnetic, and Neel-type antiferromagnetic phases, at 
appropriate coupling constant values, while on the other 
hand we find competition also with exotic spin-Peierls- 
like dimerized orders. Our results, for the SU(4) sym¬ 
metric point, are in agreement with the algebraic color 
liquid state found in Ref. [63] . 

The paper is organized as follows. In Sec. H we start 
with the generalized Hubbard model describing the four- 
component Fermi gas on an optical lattice. Then, we 
briefly summarize the steps needed to derive a superex¬ 
change model for the system in the Mott insulator limit 
with one particle per site. In Sec. HI a mean-field ap¬ 
proximation is applied to the magnetic superexchange 
model. By analyzing the mean-field solutions we charac¬ 
terize and discuss the possible ground state phases of the 
model in Sec. IV. Finally, we conclude and summarize 
our results in Sec. V. 


where the first sum runs over nearest-neighbor pairs and 
spin components a = {—3/2,—1/2,1/2,3/2}, while the 
second sum runs over the single sites of the lattice and 
spin components a and /3. We use the convention that 
Greek letters denote the z-component of the atomic hy- 
perfine spin. From now on, an implicit summation over 
repeated Greek indices is also assumed. The operators 
c\ Q, and Cj ^ respectively create and annihilate a single¬ 
particle state at site i with spin projection a. The spin- 
dependent on-site interaction is described by the tensor 
V^f ■ The total spin and its z-component of the colliding 
particles are conserved. The most general form of such 
spin dependence can be expressed with the help of pro¬ 
jection matrices that project to the two particle tensor 
product space with a given total hyperfine spin of the 
two colliding particles. 

=5o(Po)“f+ff2(P2)“:f. (3) 

Pq and P 2 are the projectors projecting to total spin 0 
and 2, respectively. These operators are antisymmetric 
under the exchange of their upper or their lower indices, 
i-e- = -iPe)jls = -iPe)s’^, where e = 0, 2. We 

note, that Pi and P 3 is missing from the sum, since they 
are symmetric in the respective spin indices and due to 
the Pauli principle their contribution cancels in Eq. (2). 
The coupling constants go and 52 are expressed in the 
usual way together with the hopping amplitude t with 
the help of the Wannier-function overlap integrals [83] . 

When the number of atoms is equal to the number of 
lattice sites, furthermore, the on-site interaction is much 
larger than the tunneling amplitude, then multiple occu¬ 
pancy becomes energetically costly. In this limit, all of 
the low energy states have exactly one particle per site. 
The dynamics restricted to this low energy sector con¬ 
tains only spin fluctuations. The superexchange Hamil¬ 
tonian, H, governing such a Mott insulator dynamics can 
be obtained with the help of perturbation theory. Here 
we follow the procedure and notations of Ref. [29] . 


II. MODEL 

We consider a system of ultracold spin-3/2 atoms on an 
optical lattice with honeycomb structure. The atom can 
be any of the alkali- or alkali-earth-metal-atomic species 
with total hyperfine spin-3/2. At low temperatures and 
for sufficiently deep optical lattices the atoms occupy the 
lowest band formed by the lowest states of the individ¬ 
ual sites, and the system is described by the generalized 



(4) 


In the following we express the projectors with the help 
of the SU(2) spin operators. For a spin-3/2 system, the 
single spin Hilbert space is 4 dimensional, and the spin 
matrices are 4x4 Hermitian matrices. We use the rep¬ 
resentation, where the single site spin basis vectors are 
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eigenvectors of Fz. The spin matrices are 


Fz = 


Fy = 


" 3 

0 

0 

0 
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0 
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2 
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n 
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3 
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73 

2 

z 
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0 -i^ 0 0 

0 -i 0 
0 i 0 
0 0 0 


The two spin tensor product space is 16 dimensional. Pq 
projects to the S = 0 (spin singlet) subspace, which is 1 
dimensional, while P 2 projects to the 5 dimensional S = 2 
(quintet) subspace. Therefore, the total antisymmetric 
sector is 6 dimensional. In the antisymmetric sector these 
two orthogonal projectors span the whole space, i.e. 

{PofJ + (^2)“f = - S^,sSp,^) = (E(-));f. 

(5a) 

Here, the notation with the superscript ’(as)’ is intro¬ 
duced as a shorthand for antisymmetrization. Further¬ 
more [29], 
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T 


(Po) 


01, (3 
7,5 




a,7'F^^(5 

= (Fi.F 2)("’^\ (5b) 


where F = {F^, Fy, Fz) is the 3 component vector of the 
spin-3/2 matrices. With a straightforward calculation 
one obtains the Hamiltonian 


(( 72 ) never dominates so strongly over the singlet one (go) 
that the SU(4) invariant interaction would become at¬ 
tractive, otherwise formation of four-particle composite 
particles would be also expected which is beyond the ap¬ 
plicability of our treatment. Accordingly, the coupling 
constant a„ is always negative. Contrary, the coupling 
ttg can both be positive or negative. Accordingly, the 
spin anisotropic interaction can be tuned from a predom¬ 
inantly ferromagnetic exchange to an antiferromagnetic 
one. 


III. MEAN-FIELD THEORY 


In order to describe the system with the help of a mean- 
field theory, we assume that both the site and bond op¬ 
erators can (but not necessarily do) have a nonzero clas¬ 
sical value. These classical values of the operators will 
be denoted by a horizontal bar above the operator. Fur¬ 
thermore, we assume that the fluctuations around these 
classical values are small: 0 « (xjj - Xi,j){Xi,j - xlj)- 

This way, Xi,jxlj ~ + Xijxlj - IXijP- We as¬ 

sume similar relations for and 

By dropping the constant terms from the Hamiltonian 
(6) and using the above approximation together with the 
definitions in Eqs. (7) we arrive to the mean-field Hamil¬ 
tonian, 


Hmi — ^rilxij cl,a^j,a F Xi,j \Xi,j\ ) 


(hi) 


0'S • Fq,^ P ’ PaP 


(hj) 


zn{n^ rij + Xi,j xlj - ni) 


tts SjS, -f-Bi ,B, 


s I ‘-’j'-’j r 



where we have introduced the following two-fermion op¬ 
erators: 


r c- 

,,a 2,0: ’ 

(7a) 

r c- 

,,o^j,o 5 

(7b) 

o,P 

(7c) 

o,P 

(7d) 


The quantities rii and describe the density and spin 
on site i, respectively. In the Mott insulator state with 
one particle per site, = I on the whole lattice. The 
operators Xi,j B^j are bond operators describing 
nearest-neighbor correlations; Xij is the SU(4) symmet¬ 
ric part, while B^^- is for correlations violating the spin 
rotation symmetry. The coupling constants in the SU(4) 
symmetric and symmetry-breaking channels are a„ = 
-t^(55o-52)/(4go52) and (go - 92 )/{^ 9092 ), re¬ 

spectively. We assume that the quintet coupling constant 


' P 01,13 ‘ F a,P |Bijj ^ 

- I). (8) 

i 

The last term can be regarded as a Lagrange multiplier 
enforcing the single particle per site constraint. The 
Hamiltonian is now quadratic in the fermion operators, 
and therefore, it can be diagonalized directly. Once the 
spectrum and eigenvectors of the Hamiltonian (8) are 
known, every physical quantity can be calculated. To 
this end, quasi-particles are introduced in terms of which 
the Hamiltonian is diagonal. Their occupation number 
is given by the Fermi-Dirac distribution function. Most 
importantly, there is a set of self-consistency equations 


= = (4,aCi.a)> 

(9a) 

Xi,3 = (4,a7.a): 

(9b) 

Si = Fo,^(cJ^^Cj^^), 

(9c) 


(9d) 


which has to be solved to find the mean fields. 

In order to solve the self-consistency equations, one 
needs to introduce a unit cell repeating periodically on 
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FIG. 1. (Color online) The choice of the unit cell. There are 
8 nonequivalent lattice points inside the unit cell, which are 
illustrated as numbers inside the circles. There are 8 differ¬ 
ent bonds running inside the unit cell and 4 bonds connect¬ 
ing sites between neighboring cells. The bonds are directed 
and are enumerated with numbers. Their orientation is rep¬ 
resented by an arrow. The coordinates in parenthesis refer 
to the position of the unit cell. The shading in the figure is 
only a guide to the eye, for locating the possible tetramerized 
clusters. 

the lattice. The choice of the unit cell is a key as¬ 
sumption in the applied mean-field calculation. In this 
system, we expect the competition between symmetry¬ 
breaking states with non-vanishing, on-site, classically 
ordered spins and states with bond-centered mean fields 
with symmetry breaking or spin-disordered nature. The 
unit cell has to be compatible with all these assumptions. 

States with classically ordered, site-centered spins does 
not require any special unit cell. In a ferromagnetic state 
all of the mean-field spins (9c) point to the same direc¬ 
tion. In a Neel-like state, the site spins are alternating 
in a checkerboard manner. Thus on the honeycomb lat¬ 
tice, even the smallest unit cell, containing two sites, can 
describe both states. The non-zero mean-field value of 
bond-centered operators signals the presence of spin cor¬ 
relations on clusters of sites. In the spin-3/2 case a com¬ 
pletely antisymmetric SU(4) singlet can be formed with 
4 sites. In the simplest case, these 4 sites occupy the 
smallest tetramer, 3 sites surrounding a central one [66]. 
In order to describe these states too, the unit cell needs 
to contain all these four sites. If we want to include also 
competition with states forming larger clusters, we need 
to introduce a further enlarged unit cell. To this end, 
we introduce a unit cell containing 8 sites, as depicted in 
Fig. I. 

We have 12 x variables and 12 x 3 B ones corresponding 
to the 12 bonds, from which 8 is connecting sites inside 
the unit cell and 4 is linking together neighboring unit 
cells. These, together with the 8x3 spin variables and 
the 8 Lagrange multipliers add up to a total of 48 com¬ 
plex and 32 real, that is overall 128 real variables. The 
solution strategy is the following. The 128 mean fields 
and Lagrange multipliers are assumed to be the variables 


we are looking for. The set of 128 equations are those 
in Eqs. (9). The right hand side of the equations de¬ 
pend on the variables implicitly through the expectation 
values of the quasiparticles and the eigenvectors of the 
mean-field Hamiltonian. In the ground state {T = 0), 
the quasiparticle occupation is 1 (0) for single-particle 
states below (above) the Fermi energy. Because of the 
jump in the occupation numbers at the Fermi energy, 
the nonlinear solver of the self-consistency equations can 
get stuck. To find the solutions we first go to finite but 
low temperatures, where the Fermi-Dirac distribution is 
more smooth. The occupation numbers are calculated 
with the help of the thermal average defined by the grand 
canonical density matrix p = exp(—/3iLmf)/Z, with the 
Hamiltonian (8), (3 = I/(fcr) the inverse temperature 
and Z = Tr[exp(—/3iLmf)] the grand canonical partition 
function. Then, we follow the solution by lowering the 
temperature to a much lower value. In our calculation the 
characteristic energy of the problem is |a„|. By choosing 
a final inverse temperature /3 = 100 x |a„|“^ we get re¬ 
sults corresponding to the zero temperature limit up to 
our numeric accuracy. 

Even though Eq. (4) is derived from the Hub¬ 
bard Hamiltonian (2), it has a much higher symmetry. 
Namely, Eq. (4) is invariant under local U(l) transfor¬ 
mations, c^-Q, —>■ which is not true for Eq. (2). 

This emergent gauge symmetry is the consequence of the 
one particle per site local constraint of the Mott insulat¬ 
ing state. The bond operators are not gauge invariant, 
they transform as, Xi,j and similarly, 

Bij —Bij-Consequently, a state of the sys¬ 
tem is characterized by the full set of the mean fields, 
{Xi,ji Si,Bij} —but in such a way, that those 
states are equivalent, whose mean fields are related to 
each other by a gauge transformation. Therefore, the 
states can only be characterized through gauge invariant 
quantities. The magnitude of the Xij bonds is gauge 
invariant. Instead of their phase, which is not gauge 
invariant, we can use the phase 4) of the Wilson loops 
H = jnje*'^, which are the products of the nonzero Xi,j 
bonds. Note, that when a Wilson loop is zero, due to the 
vanishing of a bond expectation value, one can always 
choose the bonds to be real along that loop. Inside the 
periodically repeating unit cell of our choice, there are 
four distinct elementary plaquettes, whose Wilson loops 
are: 


Hi = X 1 X 2 X 3 X 4 X 5 X 6 , (10a) 

n2 = X^X9Xl2X4X7Xll, (10b) 

ns = XeXsx/XaXmXg, (10c) 

n4 = X5 Xi2XioX2xIiX^ (lOd) 


Here the bond indices follow the convention used in 
Eig. I. However, the big Wilson loop encircling all 4 
plaquettes is inevitably real. It can be easily checked by 
taking the product of all 4 Wilson loops Eqs. (lOa)-(IOd). 
Thus, the phase of the fourth loop is just the opposite 
of the sum of the previous three. The site spins, Si, are 
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FIG. 2. (Color online) Energy per site vs coupling strength. 
On the bottom of the graph, the various phases (SB stands for 
spin bond) and approximate phase boundaries are indicated. 
The lines are guides to the eye for the coupling strength de¬ 
pendence of the energy per particle. 


gauge invariant, and therefore their expectation values 
are physical quantities. The bond spins, 'Bij, on the 
other hand, are not gauge invariant, and loop operators 
can not be introduced so straightforwardly. However, we 
can introduce a bond spin with the following definition: 



FIG. 3. (Color online) Illustration of the tt-Aux state. The 
sites are represented by red ellipses, while the homogeneous 
Xi,j bonds are drawn with blue lines. 


A. The TT-flux state 


bj 


XIj 

\Xi,j 


( 11 ) 


The bond spin is gauge invariant by construction, and 
is the last element needed to characterize the state. In 
summary, the mean-field solutions are uniquely (up to 
a gauge transformation) characterized by the absolute 
values of the 12 x bonds, the 8 site spins Eq. (9c), the 
12 bond spins Eqs. (11) and 3 of the 4 Wilson loops Eqs. 
( 10 ). 

The equivalence of mean-field states related to each 
other by gauge transformation is obvious on the physical 
state vector of the spin system, i.e. on the spin wave- 
function [1]. 


4'(Q!i,a2, ■ ■ ■,aN) = (0 




X,S.B 

mf 


( 12 ) 

The construction of the physical state vector from the 
mean fields is the so-called Gutzwiller projection. It is 
easy to see that all of the gauge equivalent states lead 
to the same physical spin wave function, apart from an 
unimportant global phase, after Gutzwiller projection [1]. 


IV. PHASE DIAGRAM 

In order to obtain the zero temperature phase dia¬ 
gram, we performed a massive search for the solutions 
of the self-consistency equations (9) numerically, starting 
from random initial configurations. The phase diagram 
of the system is determined by the properties of the low¬ 
est energy solutions. The energy of various solutions are 
plotted in Fig. 2 for different coupling constant ratios. 
We measure the energy in units of |a„| (remember that 
a„ < 0), and therefore only the ratio as/|a„| is relevant 
when studying the properties of the ground state. 


Let us start with the center of the phase diagram, 
where the couplings in the singlet and quintet channels 
are equal {go = g 2 ), and thus the spin flipping interac¬ 
tion vanishes (a^ =0). In this case, the effective Hamil¬ 
tonian (4) is SU(4) symmetric. In this high symmetry 
point, all of the low energy mean-field solutions preserve 
the SU(4) symmetry: there is neither site- nor bond spin 
order, = 0, Bij = 0. Such a state is exclusively char¬ 
acterized by the plaquette Wilson loops. For the lowest 
energy solution, we found that all of the SU(4) symmetric 
bond averages have the same magnitude, \xi,j\ ~ 0.771, 
as it is illustrated in Fig. 3. The Wilson loops are 
Hi = n 2 = Ha = n 4 « —0.21. As the Wilson loops on 
the elementary plaquettes have a uniform negative value, 
this state is a homogeneous tt-Hux state. The energy of 
this state is found to be E/N « —0.892|a„| per site. In a 
TT-flux state, time reversal symmetry is preserved, there¬ 
fore — despite the nonzero flux passing through the ele¬ 
mentary plaquettes — this state is nondegenerate, as it is 
expected from a potential universal resource state. The 



FIG. 4. (Color online) The fermion energy spectrum of the 
TT-fiux color liquid phase. 
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FIG. 5. (Color online) The temperature dependence of the 
dimensionless order parameter \xi,j\ of the yr-flux color-liquid 
phase at fls = 0. 


single fermion excitation spectrum of this homogeneous 
SU(4) state is shown in Fig. 4. Due to the complete rota¬ 
tional invariance in the 4 dimensional spin space, and also 
to the choice of our unit cell, the spectrum consists of 4 
bands, each of them are 8-fold degenerate. In our case of 
1 /4 hlling, the lowest band is completely filled. There is a 
Dirac cone touching between the fully filled lowest band 
and the empty second band. Therefore, massless Dirac- 
fermion excitations determine the low energy properties 
of the homogeneous 7r-flux state. We found that this state 
remains stable at higher temperatures, too, up to a criti¬ 
cal temperature around ksTc « 0.75|a„|. The magnitude 
of the homogeneous order parameter is plotted versus the 
temperature in Fig. 5. Close to the critical temperature, 
the order parameter \xi,j \ disappears at Tc like a square 
root of the reduced temperature, according to the mean- 
field exponent (3 = 1/2, as thermal fluctuations destroy 
the VBS order. 

We also found that this SU(4) symmetric tt-Aux state 
is robust against the presence of weak SU(4) symmetry 
violating perturbations, like a nonzero Og. It is particu¬ 
larly remarkable, that despite the Hamiltonian has only 
SU(2) global symmetry for nonzero Os, the ground state 
is still invariant under the higher SU(4) symmetry. This 
symmetry enlargement is unusual, although not unique 
in high spin systems [84]. Furthermore, we found that in 
the homogeneous tt-Aux state the order parameter Xi,j 
does not depend on the SU(4) symmetry-breaking cou¬ 
pling Os at all. As a consequence, the energy of the 
solution does not depend on Os, and from Eq. (8) its 
explicit form is given by = anJ2{t,j) IXulP ~ 

—0.771^|a„| • 12/8 « —0.892|o„|, where 12 is the number 
of the independent bonds, 8 is the number of the sites 
in a unit cell. The energy per particle of the solution is 
indicated in Fig. 2 by a solid line. 

The tt-Aux state is identical to the one reported as the 
ground state (algebraic color liquid) of the pure SU(4) 
Heisenberg model (og = 0) in Ref. [63]. Interestingly, 


in our mean-Aeld calculation the Xij parameters are ob¬ 
tained in a self-consistent way, and even without per¬ 
forming the Gutzwiller projection, the calculated energy 
of our state {E/N « —0.892|a„|) compares remarkably 
well to the one after the Gutzwiller projection performed 
with Monte-Garlo calculation: E/N « —0.894|a„| of Ref. 
[63]. Note, that the fermion spectrum. Fig. 4 also agrees 
to the one found in Ref. [63]. 

As the strength of the spin hipping interaction |as| is 
increasing, the homogeneous tt-Aux phase remains the 
lowest energy mean-Aeld solution up to the two distinct 
critical values in the ferromagnetic and Neel sides. Even 
beyond the critical value of Og, the homogeneous, SU(4) 
symmetric tt-Aux state remains a higher energy solution 
above the symmetry-breaking states, although, above 
« 0.26|a„|, the antiferromagnetic order becomes en¬ 
ergetically more favorable. On the other side, for nega¬ 
tive values of Og, aready above the critical Og’^ « —0.2|a„| 
coupling, the ferromagnetic spin exchange starts to favor 
a symmetry-breaking state. However, before the fully po¬ 
larized ferromagnetic order wins, in an extended region 
we found an intermediate state where spin and bond or¬ 
ders coexist, as it will be discussed in Section IV G. 


B. Conventional symmetry-breaking states 

Neel state: — In the Neel state both the SU(4) sym¬ 
metric bond and spin-bond expectation values are zero, 
Xij = = 0. In contrary, the site spins are nonzero. 

They have a homogeneous magnitude with maximal spin 
projection: ±3/2, but spins on neighboring sites are 
pointing to opposite directions. A nice example of the 
emerging antiferromagnetic order is shown in Fig. 6 a). 
We have seen above, that when the coupling Og becomes 
larger than the critical value « 0.26|a„|, the energy 
of the Neel state goes below that of the tt-Aux state. 
Note, that below a/.’^ iterations of the self-consistency 
equations starting from random initial mean-Aeld values 
usually converge to non-symmetry-breaking states. This 
also indicates the extreme robustness of the SU(4) pre¬ 
serving ground state even in presence of weak Og cou¬ 
plings. The energy of this classical Neel order shows the 
usual linear dependence on the exchange coupling Og: 
E^eex/N = —ag(3/2)^, as it can be observed in Fig. 2, 
too. 

Ferromagnetic state: — For large and negative ag/|a„|, 
the spins prefer parallel alignment, due to the strong fer¬ 
romagnetic coupling dominating the exchange processes. 
Accordingly, this homogeneous spin ordered state is a 
usual ferromagnetic state, in which only the spin expec¬ 
tation values Si are nonzero, and their values are inde¬ 
pendent of 1. A typical ferromagnetic state is illustrated 
in Fig. 6 b). Note, that we found the ground state to be 
fully polarized (|Si| = 3/2) as soon as the ferromagnetic 
state becomes the lowest energy state. This happens for 
«s < ~ —0.36|a„|. Similarly to the classical Neel 

state, in the fully polarized system the energy depen- 



a) 




FIG. 6. (Color online) Illustration of the Neel (a) and ferro¬ 
magnetic (b) states. The site spins are represented with blue 
arrows. 


dence on the spin exchange is E-pm/N = as(3/2)^. The 
line corresponding to this energy scaling is also plotted 
in Fig. 2. 


C. The spin-bond ferromagnetic state 

On the ferromagnetic side, between the homogeneous 
TT-flux and fully polarized ferromagnetic phases, the spin 
anisotropic and the SU(4) symmetric exchange become 
competitive. In their delicate balance, a mixed interme¬ 
diate state between and wins, where both spin 
order and valence bond expectation values on nearest 
neighbor dimers coexist. Thus, in addition to the ho¬ 
mogeneous nonzero spin averages, both of the two link 
operators Xi,j Bij have nonzero expectation values 
along the dimers. 

The above state, in which a spin anistropic dimer or¬ 
der becomes superimposed upon the ferromagnetic back¬ 
ground, has the lowest energy between and i.e. 
in the vicinity of the tt-Aux state. There are two gauge 
nonequivalent states, which are not related to each other 
by lattice symmetries either: in one case the dimers form 
a staggered columnar pattern [Fig. 7 a)], and in the other 
case, they lay in alternating directions in the neighbor- 


n ^ / 

/ iij lu it| ttf ii| 1 / 

/I ^ / 

{ li| m lu 1 

tt/ lit W lit i / 



FIG. 7. (Color online) Illustration of the columnar-bond (a) 
and zip-bond (b) ferromagnetic state. Note the alignment of 
the dimers on the two subplots. The site spins are represented 
with blue arrows, while the bond spins are shown with green 
arrows. 


ing columns forming a zip-like pattern [see Fig. 7 b)]. We 
found that this state is also robust, the magnitudes of the 
nonzero order parameters are not sensitive to the tuning 
of the spin flipping interactions: |Si| « |Bij| « \Xi,j\ ~ 
1. The mean-field energy of the solution from Eq. (8) 
is Esb/N = J2{i,j)Mxi,j\^ + -p as|Si||Sj|) « 

(an + 4as)/2, where we used that the number of dimers 
per unit cell is 4. In Fig. 2 the dependence of the en¬ 
ergy of the dimer spin-bond solution is also shown with 
a line. The energy of this dimer spin-bond state goes 
above that of the ferromagnetic state for Ug < a^'^, but 
it still remains a solution, as the magnitude of the fer¬ 
romagnetic exchange is increasing. Contrary, above the 
upper critical coupling (og > « —0.2|a„|) the SU(4) 

invariant exchange destroys spin-bond order. 

This intermediate spin-bond phase can be regarded as 
a polarized version of the SU(4)-symmetric dimer phase. 
To our knowledge, such a state has not yet been found 
in numerically exact calculations; however, in Ref. [85], 
in the framework of a slave-particle theory, a similar, po¬ 
larized quantum spin-liquid phase was found. 

In the uniform ferromagnetic part of the phase dia¬ 
gram, but not so far from the dimer spin-bond order 
(Os’^ « —0.66|a„| < Os < Og’^) we found that the compe- 
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FIG. 8. (Color online) Illustration of the disconnected-chain 
state. The Xi ,3 bonds are drawn with blue lines. Other mean 
fields are zero. 

tition between the two exchange channels is even stronger 
leading to various spin-bond states in which the nonzero 
site spin and bond averages form different complex pat¬ 
terns. The energy of these spin-bond configurations lies 
lower than that of the dimer spin-bond state. These so¬ 
lutions are not robust, and show a large diversity. 

D. Other competing states 

Above the lowest energy states, we found several other 
states, some competing with the ground state and some 
which are situated energetically well above. The compet¬ 
ing states located energetically close to the homogeneous 
TT-flux state in the < o^d region are especially 

interesting. These states are all similar to the ground 
state in many aspects: they are also SU(4) invariant VBS 
without any spin order. That is, for these states only 
the link order parameter Xi,j has nonzero expectation 
value, and they are robust against the presence of weak 
SU(4) symmetry breaking spin flipping processes. How¬ 
ever, these states are not homogeneous anymore, usually 
the translational invariance is preserved in one direction. 
Therefore these states are three-fold degenerate, corre¬ 
sponding to the lattice symmetry of rotation by 60 de¬ 
grees. We mostly found various weakly coupled ladder 
patterns with stronger bonds along the legs and rungs. 
The states with different patterns lay energetically close 
to each other. In Fig. 8 we present an extreme of such 
states, where complete dimension reduction can be ob¬ 
served: the bonds along the rungs, and between the lad¬ 
ders are zero, i.e. the nonzero bonds form disconnected 
chains. The value of the order parameter along the chains 
is \Xi,j\ ~ 0.897. As a consequence of the dimension re¬ 
duction, the fermion spectrum is flat in the reciprocal 
space along the direction orthogonal to the orientation 
of the chain. 

For tts < with energies well above the ferromag¬ 
netic, and lowest energy spin-bond states, we found sev¬ 
eral higher energy solutions of the self-consistency equa¬ 


tions where the site centered spin order coexists with 
valence bond order. In these states, the bond spins are 
disordered and only site spins exhibit the ferromagnetic 
order. 


V. SUMMARY 

In this paper we have studied the magnetic proper¬ 
ties of Mott insulators realized with ultracold fermions 
on an optical lattice with honeycomb symmetry. In the 
framework of a mean-field theory, incorporating both site 
and bond orderings on equal footing, we calculated the 
phase diagram and identified the low energy competing 
magnetic states. We found, that even though the AKLT 
parent Hamiltonian can not be realized with ultracold 
atoms, a color liquid state with a tt-Aux per plaquette 
emerges as the ground state for an extended parame¬ 
ter region surrounding the SU(4) symmetric point. This 
state is characterized by a completely homogeneous set 
of valence bonds, similar to the AKLT state. 

When the spin changing coupling constant is suffi¬ 
ciently large compared to the SU(4) symmetric coupling 
constant, the color liquid state goes to a symmetry break¬ 
ing state. In the antiferromagnetic side there is a direct 
transition to the classical Neel order from the homoge¬ 
neous TT-flux state. Contrary, in the ferromagnetic side of 
the interaction, the transition to the fully polarized ferro- 
magnetically ordered state goes through an intermediate 
phase, which is characterized both by site and bond spin 
orders. In this intermediate phase we found a narrow 
region where the emerging state can also be regarded 
as a dimerized VBS state, but with anisotropic nearest- 
neighbor correlations. For stronger ferromagnetic ex¬ 
change we found a spin-bond ’’disordered” state, in which 
the nonzero links form various disordered patterns as a 
consequence of the competition between the SU(4) sym¬ 
metry breaking(as) and preserving (a„) couplings. When 
the ferromagnetic coupling is further increased, the sys¬ 
tem ultimately goes to a ferromagnetic state with only 
site spin order. In the SU(4) symmetric point, our mean- 
field results are supported by the numerically exact meth¬ 
ods of Ref. [63], therefore we believe that at the SU(4) 
point, and also in a vicinity of it, the mean-field result 
can be trusted. Also, for high values of josj, where one 
anticipates symmetry-breaking solutions, the conclusions 
of the mean-field theory look solid. Between these differ¬ 
ent symmetry limiting cases there has to be at least one 
phase transition. As the mean-field theory is not suitable 
to describe criticality, it is hard to tell the precise loca¬ 
tion of the transition point (s). The intermediate spin- 
bond phase, between the ferromagnetic and color-liquid 
phases can not be justified, to our knowledge, by ear¬ 
lier numerically exact calculations. However, a similar 
phase was found in Ref. [85] for a different model in a 
slave-particle calculation. 

These phases can be distinguished in experiments by 
measuring quantities sensitive to the type of ordering. 
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The most striking difference between the various phases 
is exhibited in nearest-neighbor spin correlations, which 
can be measured e.g. with the superlattice technique 
[16]. Another way of testing the magnetic properties is 
by measuring the spin structure factor by spin-sensitive 
Bragg scattering [86]. An important question is how low 
the temperature should be in order to access these cor¬ 
related phases. The color-liquid phase was found to be 
robust up to a critical temperature in the order of ja„j, 
that is a few nanoKelvin in ultracold atom experiments 
on optical lattices. As for the symmetry breaking phases, 
we don’t expect them at finite temperatures in an infinite 
system due to the Mermin-Wagner theorem. However, in 
experiments in a finite system a tendency towards order¬ 


ing can happen even in finite temperatures. The zero 
temperature results can be valid up to the gap of the 
excitations, which goes approximately with jogj. 
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